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^ ■ Abstract 

> 

O ■ The Fokker-Planck equation for the probability /(r, t) to find a random walker at position r at 



time t is derived for the case that the the probability to make jumps depends nonlinearly on /(r, t). 
The result is a generalized form of the classical Fokker-Planck equation where the effects of drift, 



1-^ ' due to a violation of detailed balance, and of external fields are also considered. It is shown that 

O 

^ ■ in the absence of drift and external fields a scaling solution, describing anomalous diffusion, is only 

possible if the nonlinearity in the jump probability is of the power law type (~ f'^(r,t)), in which 
'^ I case the generalized Fokker-Planck equation reduces to the well-known Porous Media equation. 

Cd ' Monte-Carlo simulations are shown to confirm the theoretical results. 
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I. INTRODUCTION 

Random walks are typically characterized by the probability to find a walker at some 
position r at some time t, f (r, t) . This could equally well be the concentration of walkers 
at a given space-time location in the event that there is an ensemble of independent walkers. 
In the classical case of an unbiased, discrete-time random walk on a lattice, it was first 
shown by Einstein that at length and time scales large compared to the lattice spacing and 
the time step, respectively, the distribution satisfies the Fokker-Planck equation 

which is the classical diffusion equation with diffusion constant D [1]. The diffusion con- 
stant can be expressed in terms of the microscopic dynamics of the problem, namely the 
probability for the walker to take a step, the lattice spacing and the time step. It is obvious 
by inspection that the diffusion equation admits of normalized scaling solutions of the form 
/ (r, t) = t~^/^0 (^^/^) which immediately implies typical diffusive scaling of the second mo- 
ment, (r^)^ = 2Dt, which is the Einstein relation. However, there are many systems observed 
in nature where it seems natural to use the language of diffusion, but for which the mean- 
squared displacement scales as something other than linearly with time. In order to describe 
such systems, the functional form of eq. ([T]) is often generalized so as to allow for other types 
of scalings. Two popular methods are (i) the replacement of the second derivative on the 

nn 

right by a fractional derivative giving the Fractional Fokker Planck equation [2|, |3|] (FFP 
and (ii) replacing / (r, t) on the right by /° (r, t) giving the Porous Media equation [4, |5|, 
(PME). The FFPE can be understood as being the continuum limit of a continuous time 
random walk in which the waiting times obey the Mittag-LefHer distribution [7] (a general- 
ization of the Poisson distribution in which the probability for a jump decays algebraically 
with time for long times). It is therefore possible in this case to relate the mathematical 
formalism (the FFPE) to a microscopic description (power-law-distributed waiting times). 
The purpose of this paper is to describe a similar class of microscopic dynamics for which 
the PME arises naturally as the corresponding Fokker-Planck equation. 

There have been several attempts to provide some dynamical context for the PME. Abe 
and Thurner p] attempted to generalize the classical derivation of Einstein by introducing 
the concept of escort probabilities into the master equation for a random walk. Aside from 
the ad hoc nature of the generalization, the result is the PME plus an additional term which 



is not well-behaved in the long-time limit. Several authors have described the relation of 
the PME to a continuous time random walk. In particular, Curado and Nobre |9(] show 
that the PME arises from a continuous time random walk in which the transition rates, 
which are constants in the classical random walk, depend on some power of the distribution. 



Borland pi|] and Anteneodo and Tsallis [ll| discuss the fact that the PME corresponds to a 
Langevin equation with multiplicative noise but, given the equivalence of the Fokker-Planck 
and Langevin descriptions, this is just another way of writing the same result. Lutsko and 



Boon [12|] show that an assumption of nonlinear response in an ordinary fluid leads to the 
PME, but with no indication of the origin of the nonlinear response. Another approach 
based on generalizing the cumulant expansion of the intermediate scattering function leads 



to a somewhat different generalization of classical diffusion 13 1. 

Consider a discrete time random walk on a one dimensional lattice under the condition 
that the probability that the walker makes a jump from one lattice site to another depends 
on the concentration of walkers everywhere on the lattice. In this way, we generalize and 
extend previous models in several ways. First, we allow for jumps of arbitrary length and 
with asymmetric probabilities so that detailed balance is violated and an intrinsic drift is 
generated. Second, we start with a discrete-time model rather than the continuous time 
random walk which, combined with the drift, leads to additional terms in the Fokker-Planck 
equation. The continuous time random walk is a special limit of our formulation. Third, 
we do not assume a priori power-law nonlinearities as has commonly been the case. We 
allow for a quite general form of nonlinear dependence of the jump-probabilities on the 
local distributions and we show that the resulting Fokker-Planck equation only admits self- 
similar, i.e. scaling, solutions if the nonlinearities take the form of power laws. Thus we 
derive the existence of power-law nonlinearities rather than impose them. Finally, we also 
consider the effect of an external field. A preliminary discussion of these results has recently 
appeared |l4l|. 

In the next section, we start from the master equation and use a multiscale expansion 
to derive the Fokker-Planck equation. The modifications necessary to take into account 
the action of an external field are also discussed. In section III, we explore the properties 
of the generalized diffusion equation. In particular, we show that self-similar solutions 
are only possible under conditions that reduce our equation to the PME. We also present 
numerical results which demonstrate the importance of the non-standard terms occuring in 



the generalized diffusion equation and showing, in particular, the effect of breaking detailed 
balance. The last section gives our conclusions. 

II. DERIVATION OF THE GENERALIZED FOKKER-PLANCK EQUATION 
A. The master equation 

Consider a walker on a lattice whose sites are labeled by a discrete index, /. A classical 
random walk is characterized by a set of probabilities {pj} which give the likelihood for a 
jump of j lattice sites ( j > corresponds to jumps to the right, j < to jumps to the left). 
An individual walker is characterized by the probability to be at site / at time step i, fi (i). 
Equivalently, one could imagine a population of independent walkers which all start from 
the same site, in which case fi{i) would be the concentration of walkers at site I at time 
step i. If the walk is symmetric, p_j = pj, the walker exhibits diffusive behavior whereas 
asymmetric probabilities give rise to diffusion superposed on a systematic drift. 

In the present case, we generalize this picture by considering that the jump probability 
is a function of the occupation probability (or the concentration of particles) on the lattice. 
Consequently the probability to make a jump of length j from site / will depend on the 
concentration at site / at time step i, fi (i), and on the concentration at the end point of 
the jump, fij^j (i). It is convenient to introduce the more general notation whereby the 
probability to jump from site / to site k at time t is P {I —>■ k; t) so that the distribution 
obeys the master equation 

^(z + l) = ^(z)+ J2 [fk{t)P{k^l;t)-fi{t)P{l^k;t)'\ (2) 

fc=— oo 

where the first term on the right is the increase in population due to walkers jumping to site 
/ from all other sites, k, whereas the second term is the loss due to walkers leaving site / to 
go to site k. In the classical case, the jump probabilities take values drawn from a prescribed 



distribution, i.e. P {k —>■ l]t) = pi^k ISl]. Here, we make the specific generalization that the 
probabilities have the form 

P{k^l-t)=Pi_kF(jk{^Ji{€}^ (3) 

for some, as yet unspecified, function F{x,y). Note that the probabilities must satisfy the 



obvious normalization 

l = J]P(/^A:;t). (4) 

k 

This, together with the requirement that the probabihties be bounded, < P (A; ^ /; t) < 1, 
places restrictions on the form of F (x, y). 

B. Smoothing 

The goal is to examine the distribution on length and time scales that are large compared 
to the lattice spacing 6r and time step 6t. On these scales, it is expected that the distribution 
can be approximated by a continuous function. To formalize this, a smoothed version of the 
distribution is defined by 

/(r,t) = 5^G(r-Wr,t-^5t)i^W (5) 

where the sum extends over all values of the indices and the function G{r, t) is assumed to 
be localized near the point r = 0,t = 0. For example, the smoothing function could be a 
product of gaussians. 

In general, the length and time scales associated with the smoothing can be as small as 
those of the random walk model. In the following, it makes no difference as long as both 
are small compared to the scale of typical variations in the distribution function. In general, 
we assume that, as in this example, there are scales such as o"r and o"t that characterize the 
range of the smoothing and henceforth that these scales are similar to the lattice spacing 
and time step, 

l<Sr/^r,St/^t. (7) 

It will be necessary to also define the inverse transformation. To that end, notice that 
/ {kSr, j5t) = J2G {k5r - I5r, jSt - i5t) fi {i) (8) 

and assume that this relation is invertible so that 

fi{{) = Y,G-^ {k5r - I5r, j5t - t5t) f {k5r, j5t) . (9) 



This expression can be developed in a Taylor expansion. Assuming that the smoothing 
function and its inverse are even functions of their arguments, as is natural, then 



fi{i) = f {I5r, iSt) J2 G'^ (kSr - ISr, jSt - iSt) 



k,j 



+ I {5rf ^"f^^^'-f^) Y^^k- If G-' {kSr - I6r,j6t - z5t) 



2' ' (9r2 

1 ,,^,2dy{l6r,t6t) 



k,j 



+ ^ m' ' qI^ ^ E (^' - 'fG-' (^^^ - ^^'^3^' - ^^t) + - (10) 

k,j 

It is easy to see that if G {k6r — I6r,j6t — i6t) is normalized, then so is the inverse function. 
The sums then characterize the width in space and time, respectively, of the inverse smooth- 
ing functions which will be of the same order of magnitude as that of the actual smoothing 
functions. So, we have that 

/; («) = / (Wr, t6t) + 7^(T^ — + -ftcrt ^ + ... (11) 

for some dimensionless constants 7^, 7t which are of order unity. Note that this expan- 
sion makes sense, as does the whole smoothing procedure, provided the gradients of the 
distribution are small over the scales y^, \fo~t- 



C. Expansion of the master equation 

In the limit of classical diffusion, when the transition probabilities take values from a 
given distribution, one could simply multiply the master equation by G (r — /(5r, t — i6t) and 
sum to get the exact master equation for the smoothed distribution, 

CO 

f(r,t + 6t) = f{r,t)+ J2 [fir-k6r,t)p^-f{r,t)p.^] (12) 



m=— 00 



However, the nonlinearities of the generahzed model do not permit this. Instead, eg. (ITT]) is 

used to get 

d^f(r,t + 6t) dy(r,t + 6t) 
f{r,t + 6t) + 7.^.-^^^;:^ ^ + ^tat-^, ^ + ... 

oo 

+ ^ [fir- m6r, t) F{f{r- m6r, t) , / (r, t))p^ - f (r, t) F (/ (r, t) J {r - mSr, t)) p_„] 



m,=— oo 



+ lrar Y. f^^\T^'''^ Fifir-^Sr,t),f{r,t))p„ 



m=— oo 



d'f{r,t) 
g^2 



F{f{r,t),f{r-m6r,t))p_m] + 



(13) 



where we have only explicitly written one of several terms in the sum proportional to cr^ 
(and none of the terms proportional to at). The reason is that we will now further expand 
the distribution so as to give a superficially local expression. Then, it is found that the 
terms proportional to 0"^ and cxt only contribute to third order in the gradients, so that we 
have 
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df{r,t) 
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dxF {x,y) dxF{x,y) 



dx 



dy 



-1/ 
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m^Vm 



m=—oo 
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dxF (x, y) dxF (x, y) 



dx 



dy 



E 



m^Pra 



m=— oo 



d'^xF (x, y) d'^xF (x, y) 



dx"^ 



dy"^ 



f 



+ O all'' 



,d'f 



(14) 



where we have used the assumption that Sr < ^Ja\- to replace 8r by ar in the error estimate. 
A compact notation has also been introduced whereby 

' dxF {x,y)\ f dxF {x,y)\ 



dx 



(15) 



D. Multiple time scales 



We could simply truncate the expansion obtained so far on the grounds that the gra- 
dients are small over the scale of the smoothing (i.e. small over the scale of a few lattice 



spacings) but this is unsatisfactory on both physical and mathematical grounds. Physically, 
the resulting equation does not reduce to the diffusion equation in the appropriate limit of 
F{x,y) = 1. Mathematically, this results in a second order equation in time whereas the 
exact master equation is clearly first order in time: knowledge of the distribution at time 
step i is sufficient to calculate it at all future time steps. These problems are not unrelated: 
both are due to the fact that changes in the distribution in time are driven by spatial gradi- 
ents so that in some sense derivatives in time and in space are interchangeable. Ideally, we 
would like to say that the first order spatial gradients drive the first order time derivative, 
the second order gradients second order time derivatives, etc. However, we cannot simply 
equate these different terms separately as there is only one distribution and it can only 
satisfy one equation. The solution is to generalize the distribution to have many different, 
but related, time dependencies that can be satisfied at different length scales. This leads to 
the method of multiple time scales. 

To separate the different length and time scales in the problem, first define a length scale 
i over which the relative variation of the distribution is of order one, 



1 df 



or 



(9 In/ 



1. 



(16) 



/ dr/£ ' ^ " dr 

Then, a small parameter e = 6r/i is defined which quantifies the notion that the derivative 
of the distribution is small over the length-scale of the smoothing function (which we assume 
is a few lattice spacings so that 6r ~ ^/o\^)■ A parameter r is introduced by defining 6t = er 
and dimensionless variables z = r/i and s = t/r are used to write the master equation as 

^ df{z,s) ^ l^, d'f{z,s) 
ds 2 ds'^ 



df {z, s) ^ 
oz 

1 ,dy{z,s) 



+ 2^ 



^-/ 



dz^ 

dfiz,s] 
dz 



J2 



dxF (x, y) dxF {x, y) 
dx dy 

dxF (x, y) dxF (x, y) 



-1/ 



dx 



dy 



J, 



d'^xF (x, y) d'^xF (x, y) 



dx"^ 



dy"^ 



Oie^ 



(17) 



-1/ 



where J„ = X]m"^"Pni- Additional time scales are now introduced by generalizing the 
distribution to a function of many time variables f {z,s) — * f {z,So,Si, ...) where the con- 
nection between this generalized function and the actual distribution is that / {z, s) = 



/ {z, s, es, e^s, ...). Thus, the time derivatives must be replaced by 

d_ 
ds 



dsr, d dsi d 
+ 



_d_ _d_ 
dsQ dsi 



(18) 



ds dso ds dsi 

We can now demand that the terms cancel at each order in e since this just defines the 
dependence of the distribution on the various time-scales. The first two orders in e gives 



dsi 



dsQ dz 



2 dsi 



2(9^2 



J, 



dxF (x, y) dxF (x, y) 
dx dy 

dxF {x, y) dxF (x, y) 



dx 



dy 



-1/ 



-1/ 



+ 



1 fdf 



2\dz 



J2 



d'^xF (x, y) d'^xF {x, y) 



dx"^ 



dy^ 



(19) 



/ 



Now it is clear that the first equation can be used to rewrite the second derivative with 
respect to sq in terms of spatial gradients, 



9V 
dsi 



dso \ dz 



dxF (x, y) dxF (x, y) 



dx 



dy 



2 <9 fdxF{x,x)\ df 
^ dz \ dx J f 9z 



givmg 



dl 
dso 

dl 
(9s 1 



dz 



19V 

2dz^ 



J2 



dxF (x, y) dxF (x, y) 

dx dy 

dxF (x, y) dxF (x, y) 



f 



dx 



dy 



-Kl' ^^ 



d'^xF (x, y) d'^xF (x, y) 



dx'^ 



dy^ 



I j2 d f dxF (x, x^^ 
^2 ^Wz\ di /^ 



dl 
dz 



Summing gives the desired result, 

df 



ds 



dz 



dxF (x, y) dxF (x, y) 

dx dy 

dxF (x, y) dxF (x, y) 



dx 



dy 



2I9I' "^^ 



d'^xF (x, y) d'^xF (x, y) 



dx'^ 



_lj2d 
2^'dz 



dxF (x, x) 
dx 



dl 

dz 



dy"^ 
+ O (e^) 



-1/ 



(20) 



(21) 



(22) 



In terms of the original variables, this reads 



dt 



dr 



dr 



dxF (x, y) dxF (x, y) 



dx 



dy 



|/(M) 



1 



- ^C'6t 



d fdxF (x,x)\ d 



dz 



- 



dx 



dr 



f{r,t) + 0{e') 



(23) 



where 



Sr 



J2 



2 \St 
Alternatively, the diffusion coefficient can be written in terms of the second cumulant as 

1 {6r'^ 



(24) 



D = - 



2 5t 



{■h - Jl) 



(25) 



and the equation re-arranged to give 



^ + C— (xF(x,a;)) - D^ f dxF{x,y) dxF {x,y) \ 9/ 
dt dr ^ ■^ dr \ dx dy J r dr 



1 2r 9 I dxF {x,y) dxF {x,y) fdxF{x,x 

2 dr \ dx dy \ dx 




(26) 



Equations (1251 [2B]) give the generalized Fokker-Planck equation and are the main result of 
this section. 

E. Effect of an external field 

If the walkers are subject to an external field, V{r) , the derivation given above must 
be further generalized. In stochastic algorithms, such as the standard Metropolis Monte 
Carlo algorithm, the goal is to generate the canonical distribution [l5|. In fact, similar 
reasoning also lies behind the fluctuation-dissipation relation that is needed to specify the 
autocorrelations of the noise in Langevin models. Here, we can adopt the same approach 
and demand that the effect of the field be to modify the jump probabilities so as to generate 
some specified steady state distribution. Alternatively, one can adopt the position often used 
in modeling nonequilibrium processes and assume that the effect of the field is the same as 
in an equilibrium system - which would be equivalent to an assumption of local equilibrium. 
Both possibilities are explored here. 
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1. Detailed Balance 

The idea is that the stationary distribution is specified a priori as some function of the 
external field. If the stationary probability to find a walker at site fc is vr^ , then the master 
equation demands that 

oo 

= X^ [^kPi-kFki{nk,7Ci) -nipk^iFikini,nk)] (27) 

fc=— oo 

where the subscripts on the F -functions indicate that these now depend on position via the 
field. The usual condition of detailed balance would be that forward and backward jumps 
must balance, 

TlkPl-kFkl (TTfc, TTi) = TTiPk-lFik (tT/, TTfc) (28) 

However, this assumption is problematic since the elementary probabilities pi-k may make 
the forward and backward directions asymmetrical - in the extreme case, backward jumps 
might be forbidden altogether. This is simply a manifestation of the fact that asymmetric 
elementary probabilities give rise to drift and in the case of drift it makes no sense to speak 
of the stationary distribution. So we can only attempt to enforce detailed balance when the 
elementary probabilities are symmetrical, in which case fl28|) reads 

T^kFkl (vTfc, ^l) = TTlFik (ni, TTfc) (29) 

Then, making the usual separation of the jump probabilities into the probability to generate 
a particular jump, F (Tii-m, t^i) as before, and the probability to accept a jump, Gi-m,i , the 
balance condition becomes 

TTkF {TTk, ni) Gkl = TTlF (^z, TTfc) dk (30) 

which is solved, e.g., by the Metropolis ansatz 

Gkl = mm 1, ^ ^ . (31) 

V TTkF{TTk,TTl)J 

To proceed, we make the further assumption that the stationary distribution is a local 
function of the field, vr^ = $ (/?VJ) = $ {f3V (ISr)) . It is shown in Appendix lAl that in this 
case the generalized equation becomes 

% + l{c-D'{r)K {(3V (r)) ^^) {xF (x, y))^ 

2" 



d_ 
dr 



— f dxF{x,y) _ dxF{x,y) \ _ 1 ^2 e^ /' dxF{x,x) \ ' 
\ dx dy J f 2 \ dx J 



dr 



(32) 
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where 



and 



K{V) 



D'(r) 



dlnxF {x,y) dlnxF {x,y) 



dy 



dx 



HV) 



_d_ 
dV 



$(\/) 



(33) 



6t 



2 oo / f) \ 

- J2 m'pmel-mK{(3V{r))—(3V{r)y (34) 



If the elementary probabihties are symmetric, then D' (r) = D. In this case, the advection- 
diffusion equation can be written explicitly as 



a/ d_f{xF{x,y))ffdxF{x,y) dxF {x,y)\ 9$ 



dt 



Or \{xF {x,y))^ \ dx 
d fdxF{x,y) dxF {x,y)\ df 
dr \ dx dy J r dr ^ 

where the fact that / = $ is a stationary solution is obvious. 



dy 



dr 



(35) 



2. Local Equilibrium and Superstatistics 



If, on the other hand, we make the local equilibrium assumption that the acceptance 
probabilities are the same as in an equilibrium system, 

Gki = min (1, exp {-(3 {V {I5r) - V {k6r)))) (36) 

then the results in Appendix Rl give the same form as eq. fl32|) . but with K {V) = —1. 

The local equilibrium assumption can be relaxed by using the superstatistics ap- 
proach IGj better suited for systems out of equilibrium where the Boltzmann distribution 
exp {—(3 {V (r))) cannot be expected to hold. The acceptance probabilities are then written 



as 



with 



Gki = min 1, exp -(3 {U {I6r) - U (Mr)) 



rf/^^e-'^^W 



exp -pU{r] 



(37) 



(38) 



Z{f3) 
where /(/?) is a prescribed distribution of the intensive variable f3 with the normalization 

Z{P). We then obtain the generahzed advection-diffusion equation (see Appendix R|) 



dl 
dt 



+ 



C 

d_ 

dr 



d_ 

dr 



D(r) 



d (dpU{r) 
dr 



dr 



{xF{x,y)). 



^ f dxF{x,y) _ dxF{x,y) \ _ 1^2 e^ f dxF{x,x] 
\ dx dy ) f 2 \ dx 



dr 



(39) 
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with 



D(r) 



6t 



y m^PmQ m 



dpU (r) 
dr 



(40) 



III. PROPERTIES OF THE GENERALIZED DIFFUSION EQUATION 
A. The Generalized Equation as a Conservation Law 

The generahzed diffusion equation contains an exphcit velocity, C. However, since it 
multiphes a nonhnear function of the distribution, it is not a drift in the usual sense that it 
can be eliminated by a Galilean transformation. Although it arises from the same physical 
source as the drift in classical diffusion - namely, the asymmetry of the jump probabilities 
- it corresponds in the present case to a position-dependent velocity. Of course, one could 
always make a transformation to an arbitrary moving frame, say with velocity C", and this 
would introduce the usual term C'drf into the equation. 

The generalized Fokker-Planck equation can also be cast in the usual form of a conser- 
vation law. 



dt 



d_ 
dr 



J = 0, 



(41) 



with flux 



J 



C-D' (r) K {(3V (r)) ^^^i ) {xF {x, y)) ^ 



dr 



— / dxF {x,y) dxF {x,y) 
\ dx dy 



_ \(j2^^ f dxF{x,x] 
2 \ dx 



dr 



(42) 



Indeed, one interpretation of the result is that it describes ordinary diffusion with drift 
velocity C and diffusion constant D that are functions of the distribution, i.e. 



^ + ^(c{r)-D'{r)F{fJ)K{f3V{r)) 



d(3V (r) 
dr 



d df 

/=|^D(r)^, (43) 

dr dr 



where 



C = CFix,y) 

' dxF (a;, y) dxF (x, y) 



D = D 



dx 



1^2^^ ( dxF{x,x) 
2 \ dx 



dy , f ^ V -- // 

This makes clear that in the special case F{x,y) = 1, classical diffusion is recovered. 



(44) 
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B. Scaling solutions 

We now specialize to the case that there is no drift, C = 0, and no external field, and ask 
under what circumstances a scaling solution of the form / (r, t) = t~"'/'^(j) {r jV'l'^^ is possible: 
in other words, when does the general formulation describe diffusion? Without drift and 
without external force, the generalized diffusion equation reduces to 

where we have introduced M(/) = ( ^ q^''^' ^ ^ '^^ j . Defining C, = r/P^"^ and introduc- 
ing the scaling ansatz gives 

-'{{hQ^C^Hq) -^'■-^(^'('-'^(0)^^(0) . (46) 

It is only possible to eliminate the factors of the time if M (/) = rriof^, for some constant 
mo, giving 

_2i^C^(C) = ™.15t'--/=A(^,(Oi^^(C)). (47) 



So scaling works provided 



7 = , (48) 



and the equation for the scaling function is 



or 



"»4(^'<^'5^^'^')^i^^^<^' = °' '''^ 



m<,D4," (C) ^'^ (C) + lC<l> (!■) = A (50) 



for some constant, A. In the case A = 0, the particular solution is easily found from 

dC^ y^) 2moD 



^0'' (C) = -^z?7n C> and is given by 



(J){0=(b ^ =C'^ '&(b ^ =C=^V (51) 

The constant B is determined by normalization: 

V 8moD ; v^ 2; ^ ' 
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where B{x,y) = T {x)T (y) /T {x + y) is the beta function. The distribution can also be 
written as a g-Gaussian by defining rj = 1 — q, 

1-^ A2^^oA 1-^ a2 



0(C) = Bi^ 1 :=^ C] 1 -=^ C • (53) 

V 2BmoD{3-q) J \ 2BmoD {3 - q) J 

Note that the scahng hypothesis demands that 

M^f^^(?5E^^^(E2yl'j =„„/., (54) 

and the fact that the function F is defined in terms of the jump probabihties means that it 
must be bounded. From its definition, we expect that 

0<xF {x, y) <l and < yF {x, y) < I , (55) 

for all x,y E [0, 1]. If for example F{x, y) = F{x), then the scaling hypothesis is: F{x) ~ x^, 
so that the bounds given above demand that 

r/>0, 0<7<1, and g<l. (56) 

All of the preceding concerning the scaling behavior only applies to the case that the 
constant A is taken to be zero in eq. fl50l) . For values oi A y^ 0, no general solution of this 
equation could be found. However, note that if {() is analytic at C = 0, then from eq. (l50|) 



lim-^0(C)= -^ , , (57) 

so that any solution with A 7^ is not symmetric about C = 0- Thus, we can say that the 
scaling behavior discussed here applies to the general case of symmetric solutions. 



IV. NUMERICAL TESTS 

In order to test the validity of the generalized diffusion equation, we have performed 
numerical simulations of the underlying random walk model. Our simulations begin with 
a population of A^, independent random walkers at position r = at time t = 0. At time 
step i, each walker makes a jump of m lattice sites from its present position, say site /, with 
a probability: pmF I fi (i) , fi+m (i) ) where the distribution fk (i) is simply the fraction of 
walkers at site k at time step i. All of the simulations discussed below were performed using 
a population of size A^ = 10^. 
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FIG. 1: (Color online) The evolution of an initial delta-function distribution for the case q = 0.999 
and equal elementary probabilities for jumps up to length 2. The symbols are from Monte Carlo 
simulation of the random-walk and the solid lines are the analytic g-Gaussian solution (j58p to the 
generalized diffusion equation. 

The first simulation is for the case of no drift, elementary probabilities pj = | with 

j G [—2, 2], and -F(x, y) = x^^'^ for which the theory gives 

/(r,t) = t-3^0g (^r/t3^) 



MQ 



B„ 






1-g . 

~ 2B^{2-q){3-q)D^ 

g)(3-g)\^ „ / 1 1 



1-9 



e 1- 



l-q 



^C 



2B^ (2 -q){3-q)D 



B 



2-2q 



(58) 



8{2-q)D J \l-q 2^ 

As stated above, only the range g < 1 is permitted, and the value q = I corresponds to 
classical diffusion. Since the initial condition and jump probabilities are symmetric, there is 
no drift and the scaling solution applies. Figures ([I]l3]) show the analytic results, given by 
eq. (l58ll . and the results of the microscopic simulations for q = 0.999 (essentially the classical 
case), q = 0.5 and q = 0.0. These correspond to anomalous diffusion with scaling exponent 
7 = 0.999 5, I and | respectively. In all cases, the agreement between the simulations and 
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FIG. 2: (Color online) Same as Fig. [H but for q = 0.5. 

the scaling solution is very good, even at the earliest times. 

In the second set of simulations, particles are subjected to drift. In this case, the elemen- 
tary probabihties are taken to be pj = ■^ for j E [—2, 2]. Figures dUE]) show the evolution 
of the distributions for the same values of q as for the no-drift case. As mentioned above, 
it is then no longer possible to solve the generalized diffusion equation analytically. So 
comparison is made to a numerical solution of eq. (!26|) with F{x,y) = x^^"^. The numerical 
solution was performed using centered finite differences in the spatial variable and a simple, 
first-order scheme in the time, with the lattice spacing fixed at 5r and the time step equal 
to 0.001(5t. For q = 0.999, the process is essentially that of the classical case of advection- 
diffusion. For the smaller values of q however, the distribution is very different, becoming 
increasingly asymmetrical as time progresses. As q becomes smaller, and the processes be- 
comes more sub-diffusive, the velocity of the peak of the distribution also decreases. Even 
with these significant qualitative changes for decreasing values of q, the generalized diffusion 
equation is again seen to give very good agreement with the Monte Carlo simulations. 

One interesting question is whether the new terms appearing in the generalized diffusion 
equation ( l23l 1261) play any role, or whether they could be neglected giving a result closer to 
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FIG. 3: (Color online) Same as Fig. [H but for q = 0.0. 



the Porous Media Equation \4^ which (with drift term) reads 



(59) 



To investigate this, we repeated the solution of two modifications of the generahzed diffusion 
equation. In the first case, the "extra" terms are simply omitted from eq. fl23|) giving 

'dxF {x,y) dxF {x,y)\ df 



dl 
dt 



d_ 
dr 



+ C^ixF{x,y))f 



dr 



dx 



dr 



(60) 



dy yf 

One objection to this approximation is that it does not reduce to the expected result in the 
limit of classical diffusion, since then the extra term would combine with the diffusive term 
to make the replacement D ^ D. This leads to the second modification considered here, 
namely omitting the extra term from equation (!26l) which then reads 



df d d f 



d fdxF{x,y) dxF {x,y) 



dx 



dl 
dr 



(61) 



dy yf 

For want of better terms, these will be referred to as modifications I and II respectively. 
Figures (JTj) and (IHl) show the numerical solution of these equations compared to the simula- 
tion data for the two cases q = 0.999 and g = 0. For q = 0.999, the type II modification is a 
much better approximation to the data than is the type I modification as might be expected 
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FIG. 4: (Color online) The evolution of an initial delta-function distribution for the case q = 0.999 
and pj = ^^jr- for j G [—2,2]. Since the probabilities violate detailed balance, there is a non-zero 
drift velocity, C = -^. The symbols are from Monte Carlo simulation of the random- walk and the 
solid lines are the numeric solution of the generalized diffusion equation, eq. (f26]l . 

since type II becomes exact for q = 1. However, the results for g = are exactly reversed: 
type I is a noticeably better approximation than is type II. The conclusion is that the full 
equation is necessary to provide a good description of the system for all values of q. 



V. CONCLUSIONS 

We have shown that a simple modification of the classical random walk gives rise to 
sub-diffusive behavior. The required modification is that the probability to make a jump 
from one lattice site to another depends on the occupation probability of the walker on the 
lattice. Using a multiscale expansion of the exact master equation, we derived a generalized 
Fokker-Planck equation. In the limit of symmetric probabilities to jump left and right, this 
equation gives rise to diffusive behavior of the moments of the distribution provided the 
dependence of the jump probabilities takes the form of a power law. In this case, our result 
reduces to the well known Porous Media equation. Unlike other approaches that begin 
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FIG. 5: (Color online) Same as Fig. U but for q = 0.5. 
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FIG. 6: (Color online) Same as Fig. [H but for q = 0.0. 
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FIG. 7: (Color online) Simulation data for q = 0.999 and the predictions of the Fokker-Planck 
equation with the type I modification. The type II is not shown as it gives virtually the same 
result as the full Fokker-Planck equation, as shown in Fig. UJ and is in almost perfect agreement 
with the data. 

with a continuous time random walk, we specifically consider a microscopic model in the 
hydrodynamic limit of large length and time scales. This is responsible for the appearance of 
a new term in the generalized diffusion equation which, as comparison to simulations of the 
microscopic model shows, is necessary to correctly describe the evolution of the distribution. 
Our generalized equation reduces to previous results in the appropriate limits. Most sim- 
ply, if the function F{x,y) = 1 the Fokker-Planck equation becomes the classical advective- 
diffusive equation. The continuous time random walk results from the scaling 



6t 



e^6t 



Sr — > eSr , Ji — > eJi 



(62) 



and the limit e — *> 0. (Note that this limit is easily deduced directly from the smoothed mas- 
ter equation without need for the multiscale expansion.) With the further approximations 
of (i) no drift ( C = 0) and (ii) hops of only one lattice site our result agrees with that of 
Curado and Nobre |9|]. 

We have shown that exact self-similar solutions of the generalized diffusion equation 
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FIG. 8: (Color online) The simulation data for q = together with the predictions of the Fokker- 
Planck equation with the type I (full line) and type II (broken line) modifications. 

(without drift) are only possible if the jump probabilities scale as power laws. In this 
case, the distribution turns out to be the so-called g-Gaussian often introduced in an ad 
hoc manner to describe anomolous diffusion. The model presented here therefore gives 
one answer to the question "what underlying dynamics could give rise to the observed q- 
Gaussian distributions": a dependence of the jump probabilities on the local distribution (or 
more hkely, local concentration) of walkers is sufficient. Note that the dependence need not 
be an exact power-law: it is enough that the long-time limit of the diffusion equation admits 
scaling which in turn implies that the function F (x, y) become algebraic in x in the limit 
of very small or very large x depending on the various scaling exponents. One restriction 
of the exact scaling result, however, is that our model is only well defined if F{x,y) = x^ 
for 7] > which in turn implies sub-diffusive scaling of the moments. To describe super- 
diffusion there are only two possibilities. Either one could construct a function F (x, y) that 
gives the proper normalization of the jump probabilities and that gives super-diffusion in the 
long-time limit or one could modify the basic description of the jump probabilities, eq.(l3]) 
so as to introduce nonlinearity in some other way. 

The generalized diffusion equation (GDE) is in some ways similar to the fractional Fokker- 
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Planck equation (FFPE): both describe sub-diffusion and both require power-law probabil- 
ities to give the subdiffusion (the GDE in the jump probabilities, the FFPE in the waiting 
times). It is natural to ask whether, given some experimental data which shows subdiffu- 
sion, there is any way to choose between the two descriptions. On physical grounds, the idea 
of waiting times that are distributed as a power law might be more appropriate in which 
case the FFPE should be prefered; if there it makes more sense to think in terms of an 
interaction between the walkers, then the GDE might be more appropriate. Empirically, if 
the distribution of walkers is measured, it might be possible to choose a model based on the 
fact that the GDE predicts that the distribution of walkers in a system showing subdiffusion 
with no external forces should obey a g-Gaussian distribution whereas in the case of the 
FFPE there is also a scaling solution but the distribution is a stretched gaussian[2]. In fact, 



in two studies, one of sub-diffusion induced by a random walk on a Sierpinski gasket 17l | 



and the other of super-diffuson induced by a raondom walk on a tree structure p^, it was 
shown that the FFPE and GDE results were sufficiently different as to allow an empirical 
distinction to be made. 

In the presence of either drift (ie non-symmetric jump probabilities) or an external field, 
the GDE is more complex than the equilvalent extension of the Porous Media equation. 
This is true even when a power-law dependence of the jump probabilities is assumed since 
in this case, the GDE becomes 



'f- ^^c+(i+„)/>'w "°^r'^" V"" 



dt dr \^ ' ^ ' '^ ^ ^ dr 



d_ 
dr 



{l + r^)DP-]^{l + ilfCHtf'^ 



dl 
dr 



(63) 



where 



r./.N (-^O ^ 2 r.f.. ^ d\n^ {13V {r))\ 

D'{r) = ^ J2 mVme((l + r/)m ^£-^21] (64) 



or, with 1 + rj = a, 

|4(c + a.Z>'(r)£MfKW)),«.^ 

ot or \ or J or 



D + |C^ (1 - ar-') 



^-r , (65) 



(9r' 

to be compared with the PME, eq.( l59l) . With no field, the drift does not generate a simple 
Galilean transformation of the equation without drift, as is usually assumed to be the case 
with the Porous Media equation, but instead generates new nonlinearities in the GDE. 
Because the drift term has the same number of powers of / but one fewer derivative, than 
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the right hand side, no simple scahng solution is evident. In the case of an external field 
but no drift, one has that D' (r) ^ D so that the gradient terms on the left and right 
hand sides of the equation have the same number of powers of / and of spatial gradients. A 
scaling solution would then be possible, but only with a trivial external field. This superficial 
analysis suggests that exact scaling is only possible in the GDE in the case of no field and 
no drift. It leaves open the possibility of approximate scaling in the long-time limit, not 
to mention the possibility that more complex assumptions for the dependence of the jump- 
probabilities might give completely different scaling properties. These questions are the 
subject of on-going research. 
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APPENDIX A: ROLE OF THE EXTERNAL FIELD 

The master equation is 

oo 

fir,t + St)=f{r,t)+ J2 Pr 

m=—oo 

where 



f {r — m6r, t) F {f {r — m6r, t) , / (r, t)) Gr- 
-f (r, t) F (/ (r, t) , / (r + m6r, t)) Gr,r+. 



(Al) 



' ™'' <I>,„„F ($,_„, $0 V <^l-rnF{<!>l-m,<^l)J \<!>l-mF {<!>l-m, <^l) J ^ ' 

and 

$; = $(V(Wr)). (A3) 

This can also be written as 



^<I);„„F ($;„„,$; 
1 + Hi^rn,l 
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where 



with 



H, 



I— ml 



h {I6r, (/ — m) 6r) 
h ((/ — m) 6r, I6r) 



h{{l — m) 6r, I6r) 



h (x, y) = <^{V (x)) F ($ {V (x)) , $ {V (y))) . 



(A5) 



The goal is to develop the expansion of Hi^rn,i in terms of 6r and to use this to derive the 
modified advection-diffusion equation. In the Appendix, we use an abbreviated notation 



whereby dr 



dr ' 



etc. 



First, note that for present purposes we need the expansion of Hi^„i,i up to order {6r) 
inclusive. For the step-function part, we have 

h{x,x — u) 



e 1 



h{x — u,x) 



& {h{x — u,x) — h {x, X — u)) 



XX ~ ^yy) 



Q{u{hy- h^) + -u^ {h, 

( u 1 

( -r-, (/''y - /^-x) + - \u\ {Kx - hyy) + 



(A6) 

Now, we assume that {hy — h^) is of order one, so that in some formal sense we can expand 
to get 



9 1 



h{x — u,x)J \\u\ / 2 \\u\ / 

In general, the (5-function (and higher order terms) only contribute on a set of measure zero 
and can be neglected. (Furthermore, we will explicitly show that the delta function cannot 
contribute until at least order {6r) .) Expanding the coefficient of the step function gives 



h{x,x — u) 



1 



u 



IIX lly 



1 



h 



_uU -^-2 



h^h 



x'^y 



2^ ^- 



+ ... 



(A8) 



h{x — u,x) J h 2 \ h K^ ' K^ h 

Multiplying these two contributions, we see that the (5-function first appears at order m^, 
but in the form of x5 {x) which is always zero; so, as stated above, it cannot contribute until 
at least order -u^, if at all. The result is 



h{x,x — u) -.\ r>, f -I h{x,x — u) 



h{x — u,x) 



h{x — u,x) 



h ) 2 \ h K^ h^ h 



e{u{hy-h^)) + 0{u^) , (A9) 
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and consequently 

Hl-m,l ■ 



{m6r) ( ^^ ) + i imSrf ( ^ " 2^ + 2^ " ^ 



Q{m{hy - h^)) . 
(AlO) 



Similarly 



Hi,i+^ = ( ^'J:z^l^T' ^!^ -1)0 (^'^ ("^^ ^^i+m) - ^i+^F ($,+„, $0) 



/i ((/ + m) 6r, I6r) 



10 1 



h{{l + m) 6r, I6r) 



h {I6r, {I + m) 6r) J \ h{l6r,{l + m) 6r) 
= (Am5r) ^^^^ + i im5rf (^ - ^ + ^ " x)) ® ^^"^^'^ ^^' " ^'^^ 
+ O (m^) (All) 



Substituting back into the master equation gives 
f{r,t + 5t) = f{r,t) 

oo 

+ J2 Pmifir- rn6r, t) F {f {r - m6r, t) , f (r, t)) - f (r, t) F (/ (r, t) , / (r + m6r, t))] 

m=—oo 

oo 

+ X^ Pm[fir-m5r,t)F{f{r-mSr,t)J{r,t))Hr-rn,r 

- f (r, t) F (/ (r, t) , / (r + m5r, t)) //,,,+„] . (A12) 



m=—oo 
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The last term on the r.h.s. is 

oo 

/(r,t)F(/(r,t),/(r,t)) ^ p„ [if,_„,, - if, 



r,r+mj 



+ 5r{drf) J2 



mp„ 



m=— oo 

dxF 
-H, 



dx 



r—m,r 



dxF 
dy 



-H, 



r,r+m 



+ ... 



/(r,t)F(/(r,t),/(r,t))x 



X 



oo _. 



^yy o hxhy , cyh^ hxx 

h ^~fP~ ^ ^P" h 

■X: 

h 



hxx "ihxhy I "ihy hyy 

h"^ ^ h"^ h 



+ {Srf{drf) Y. 



rn^Pr, 



m=—oo 



L/'eAy J^ f f v ^' i Ijq I I i-AjtA-J J~ / i V /y i Ij't i 



dx 



h 



dy 



h 



Q{m{hy-h^)) 

Q{rn{fy-U)) 



+ {5rf (dj) {^^^ 



dxF dxF 
dx dy 



^-^ + ^-^) E m'Pme{m{fy-f.)) 

^ m=—oo 
oo 

oo 

Y m'^Pni'S) (m {fy - /^)) 

i=— oo 
oo 

{6rf Y 'rn^PmQ ijn {fy - f^)) 



^^ f{r,t)F{f{r,t)J{r,t))^^ 



m=— oo 
(fhih Shih 



{drf) 



dxF _|_ dxF 
dx dy 



{5r 



d\nh dlah 

dy dx 

dlnh d\nh 



Since this term only contributes to the master equation at order {6r) , it is easy to see that 
the complete Fokker-Planck equation now reads 



dtf + Cdr {yF {y))f = D' (r) dr \ f (r, t) F (/ (r, t) , / (r, t)) 



dlnh dlnh 



dy dx 



dF{yy 
dy /^ 



1 

— ( 
2 



dF\ 



+ Dd.[F{y)-y:^^^] dJ-^6tC%{F{y)-y^+^^y^ 



where 



{5r) 



2 oo 



D'{r)=^—^ J2 m^prr,e{m{hy-h^)) . 



x=y=r 
2^ 

drf , 
f 

(A13) 



(AM) 



m=— oo 
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Note that in the case of symmetric elementary probabihties 



St 



2 oo 



^ m'^PmQ (m {fy - /^)) 



bt 



2 oo 



^ m^ {pmQ (m (/j^ - /^)) + p_„0 (-m (/i^ - h^))) 



6t 



m>0 

2 oo 



y^ "^^Pm (0 ("^ ify - fx)) + {-m {hy - h^))) 



m>0 






m>0 



2 5t ^-^ 



(A15) 



m,=— oo 



The final form of the advection-diffusion equation can be clarified. Writing it as 



dr 1 h{r,r) \dy dx J ^^y^^ 



dy 



f 



and noting that 



>c=. .,„-.^.(^ 



—h {x, y) = — $ {V (x)) F ($ {V (x)) , $ {V {y))) 
ox ox 

9$ f dxF {x,y) 
dx V dx 



drf : 



(A16) 



(A17) 



gives 



dtf + CdriyFiy)), = D'ir) 



d_ 

dr 



fFjfJ) f dxF{x,y) _ dxF{x,y) 
$F(<I>,$) V dy dx 



a$ 



+ dr 



D 



(.(.,-.^)^-fc=a„(..,-.^.(MM)'^)^ 



The local equilibrium result can also be easily deduced. It corresponds to taking 

d dV (x) 

h{x,y) = exp{-pV{x)), i.e. h{x,y) = -P ^^exp {-PV (x)) . 

ox ox 



drf. 

(A18) 
(A19) 



Noting that the derivative of D'{r) produces a term of the form x 6 (x) which vanishes, D'{r) 
can be taken under the derivative dr in flAlSp . and in the local equilibrium case 



^' ('■) <! ("^^^ (^- ^') 4 («' w ^'-^f^ (^- ^) 



(A20) 
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The resulting generalized Fokker-Planck equation reads 
dtf + d, 

dF{y) 



3V (r 
C{yF{y))^-D'{r)(p^^fF{fJ) 



+ dr 



D{F{y)-y- 



dy 



dr 
5t 



-C'\F{y)~y 



dF{y) 
dy 



f dyFjy) 
V dy 




drf. 



(A21) 
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